Home Range Crash Course

---

1 Quick Info

Coyotes (Canis lantrans) data on three individuals near Albany, NY, during the spring and summer of 2001, were used as examples for performing home range analyses.

Data comes from:

Bogan, D.A. 2004 Eastern coyote (Canis latrans) home range, habitat selection, and survival rates in the suburban Albany Pine Bush landscape of New York. MA Thesis. State University of New York at Albany.

2 Dataset

coyote<-read.csv("data/coyote.csv")
head(coyote)

Convert Lat & Long to UTM Because I’m stubborn

zone18<-data.frame(coyote$location.long, coyote$location.lat)
names(zone18)<-c("X", "Y")
zone18<-as.matrix(zone18)
UTM<-project(zone18, "+proj=utm +zone=18 ellps=WGS84")
UTM<-data.frame(UTM)
names(UTM)<-c("utm.easting", "utm.northing")
coyote<-cbind(coyote, UTM)

Lets see if the dataset has any outliers, and get a rough view of the locations in a map.

map1 <- ggplot() + geom_point(data=coyote, 
                                   aes(utm.easting, utm.northing,
                                       color=individual.local.identifier)) +
                        labs(x="Easting", y= "Northing") +
                        guides(color=guide_legend("Identifier"))

ggplotly(map1)

NA

Here, I will do a quick map on the three coyote individuals with a basemap from Open Street Map.

utm_points <- cbind(coyote$utm.easting, coyote$utm.northing)
utm_locations <- SpatialPoints(utm_points, 
                 proj4string=CRS("+proj=utm +zone=18 +datum=WGS84"))
proj_lat.lon <- as.data.frame(spTransform(
                utm_locations, CRS("+proj=longlat +datum=WGS84")))
colnames(proj_lat.lon) <- c("x","y")
raster <- openmap(c(max(proj_lat.lon$y)+0.01, min(proj_lat.lon$x)-0.01), 
                  c(min(proj_lat.lon$y)-0.01, max(proj_lat.lon$x)+0.01), 
                  type = "bing")
raster_utm <- openproj(raster, 
              projection = "+proj=utm +zone=18 +ellps=WGS84 +units=m +no_defs")


autoplot(raster_utm, expand = TRUE) + theme_bw() +
  theme(legend.position="bottom") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_point(data=coyote, aes(utm.easting,utm.northing,
             color=individual.local.identifier), size = 3, alpha = 0.8) +
  theme(axis.title = element_text(face="bold")) + labs(x="Easting",
        y="Northing") + guides(color=guide_legend("Identifier"))

Both Beta and Omega had overlaping locations.

files <- list.files(path = ".", pattern ="*.csv" , full.names = TRUE)

3 Home Range Analysis

Three basic types of home range analyses we will perform in this exercise: Minimum Convex Polygon (MCP), Kernel-Density Estimation (KDE), and Brownian Bridge Movement Model (BB).

MCP - draws the smallest polygon around points with all interior angles less than 180 degrees. MCPs are common estimators of home range, but can potentially include area not used by the animal and overestimate the home range

KDE - calculates the density of features in a neighborhood around those features.

BB - models movement of an individual

3.1 MCP Analysis

mcp_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm.easting)
  y <- as.data.frame(data$utm.northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=18 +ellps=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  mcp.out <- mcp(xy, percent=100, unout="ha")
  mcp.points <- cbind((data.frame(xy)),data$individual.local.identifier)
  colnames(mcp.points) <- c("x","y", "identifier")
  mcp.poly <- fortify(mcp.out, region = "id")
  units <- grid.text(paste(round(mcp.out@data$area,2),"ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  mcp.plot <- autoplot(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=mcp.poly, aes(x=mcp.poly$long, y=mcp.poly$lat), alpha=0.8) +
    geom_point(data=mcp.points, aes(x=x, y=y)) + 
    labs(x="Easting (m)", y="Northing (m)", title=mcp.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  mcp.plot
}

pblapply(files, mcp_raster)

  |                                                  | 0 % ~calculating  
  |=================                                 | 33% ~01s          
  |==================================                | 67% ~00s          
  |==================================================| 100% elapsed=01s  
[[1]]

[[2]]

[[3]]

3.2 Kernel-Density Estimation

kde_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm.easting)
  y <- as.data.frame(data$utm.northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=15 +south +ellps=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=100)
  ver95 <- getverticeshr(kde, 95)
  ver75 <- getverticeshr(kde, 75)
  ver50 <- getverticeshr(kde, 50)
  kde.points <- cbind((data.frame(data.proj@coords)),data$individual.local.identifier)
  colnames(kde.points) <- c("x","y","identifier")
  kde.poly95 <- fortify(ver95, region = "id")
  kde.poly75 <- fortify(ver75, region = "id")
  kde.poly50 <- fortify(ver50, region = "id")
  units <- grid.text(paste(round(ver95$area,2)," ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  kde.plot <- autoplot(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=kde.poly95, aes(x=kde.poly95$long, y=kde.poly95$lat), alpha = 0.4, fill="red") +
    geom_polygon(data=kde.poly75, aes(x=kde.poly75$long, y=kde.poly75$lat), alpha = 0.4, fill="purple") +
    geom_polygon(data=kde.poly50, aes(x=kde.poly50$long, y=kde.poly50$lat), alpha = 0.4, fill="darkblue") +
    geom_point(data=kde.points, aes(x=x, y=y)) +
    labs(x="Easting (m)", y="Northing (m)", title=kde.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  kde.plot
}
pblapply(files, kde_raster)

  |                                                  | 0 % ~calculating  
  |=================                                 | 33% ~01s          
  |==================================                | 67% ~01s          
  |==================================================| 100% elapsed=01s  
[[1]]

[[2]]

[[3]]

3.2.1 Brownian Bridge Movement

A <- read.csv("Alpha .csv")
date <- as.POSIXct(strptime(as.character(A$timestamp),"%Y-%m-%d %H:%M:%S", tz="Asia/Bangkok"))
A$date <- date
A.reloc <- cbind.data.frame(A$utm.easting, A$utm.northing,
                                as.vector(A$individual.local.identifier),
                                as.POSIXct(date))
colnames(A.reloc) <- c("x","y","id","date")
trajectory <- as.ltraj(A.reloc, date=date, id="Alpha")


sig1 <- liker(trajectory, sig2 = 58, rangesig1 = c(0, 5), plotit = FALSE)
opha.traj <- kernelbb(trajectory, sig1 = .1351, sig2 = 58, grid = 95)
bb_ver <- getverticeshr(opha.traj, 95)
bb_poly <- fortify(bb_ver, region = "id", 
                   proj4string = CRS("+proj=utm +zone=15 +south
                                     +ellps=WGS84 +units=m +no_defs"))
colnames(bb_poly) <- c("x","y","order","hole","piece","id","group")

# bb_image <- crop(opha.traj, bb_ver, 
#                  proj4string = CRS("+proj=utm +zone=15 +south 
#                                    +ellps=WGS84 +units=m +no_defs")) Crop Does not want to work

bb_units <- grid.text(paste(round(bb_ver$area,2)," ha"), x=0.85,  y=0.95,
                      gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
bb.plot <- autoplot(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
  coord_fixed(xlim = c(min(A$utm.easting), max(A$utm.easting)), ylim = c(min(A$utm.northing), max(A$utm.northing)))+
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_tile(data=opha.traj, 
            aes(x=opha.traj@coords[,1], y=opha.traj@coords[,2],
                fill = opha.traj@data$ud)) +
  geom_polygon(data=bb_poly, aes(x=x, y=y, group = group), color = "white", fill = NA) +
  scale_fill_viridis_c(option = "inferno") + annotation_custom(bb_units) +
  labs(x="Easting (m)", y="Northing (m)", title="Alpha") +
  theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5))
bb.plot

3.3 Animate Trajectory Data

Here is my attempt at creating an animation of the individual “Alpha”’s movement from April to July in 2001 in Albany, NY, USA. Data was collected by radio telemetry.

animate_frames(frames, fps = 5, overwrite = TRUE,
               out_file = "./moveVis-5fps.gif")